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Abstract. We study wave propagation in periodic and frequency dependent materials. The 
approach in this paper leads to spectral analysis of a quadratic operator pencil where the spectral 
parameter relates to the quasimomentum and the frequency is a parameter. We show that the 
underlying operator has a discrete spectrum, where the eigenvalues are symmetrically placed with 
respect to the real and imaginary axis. Moreover, we discretize the operator pencil with finite 
elements and use a Krylov space method to compute eigenvalues of the resulting large sparse matrix 
pencil. 

1. Introduction. In this paper, we study wave propagation in a two-dimensional 
periodic medium, where the material components depend on the frequency. The 
underlying operator has a band-structure, which by appropriate choice of the material 
configuration has gaps [H| . Such structures have numerous applications in optics, 
photonics and microwave engineering [14] . Electromagnetic wave propagation is in 
general frequency dependent, since dispersive effects such as resonances or relaxation 
processes are present in the material [31[T3]. The frequency dependence of the material 
is in most cases disregarded in spectral analysis. This is justified in some cases where 
the frequency dependence is weak. However, the frequency dispersion can in general 
not be ignored. 

In the case of frequency independent materials, considerable mathematical prog- 
ress has been made (24j [17] • In the frequency dependent case, the spectral problem 
is non-linear, which complicates the analysis. There exist two different approaches to 
handle the frequency dependence. Within the first method, a real-valued quasimo- 
mentum (the Floquet-Bloch wave vector) is specified and the spectral parameter is 
related to the time frequency [30, 28J . The second option is to study the quasimomen- 
tum A; as a function of a real frequency us. The non-linearity is arbitrary within the 
first approach, but the second approach leads regardless of the frequency dependence 
in the material model to a quadratic non-linearity. The quadratic non-linearity sim- 
plifies the spectral analysis and the numerical calculations. However, only numerical 
simulations exist for this approach [571 E] • 

We use the second option to study wave propagation in a periodic and frequency 
dependent medium that for a frequency interval is characterized by a space- and 
frequency-dependent real- valued permittivity e(x,us). The waves propagate in a non- 
magnetic material with the permittivity e(x\, X2, us) independent of the third coordi- 
nate x%. Let E and H denote the electric and magnetic waves, respectively. The X3- 
independent electromagnetic wave (E, H) is decomposed into transverse electric (TE) 
polarized waves (E\, E2, 0, 0, 0, H3) and transverse magnetic (TM) polarized waves 
(0, 0, E 3 , Hi, H 2 , 0) [TS]. This decomposition reduces the spectral problem for the 
Maxwell operator to one scalar equation for H3 and one scalar equation for £3. We 
consider TM polarized waves. This gives a Helmholtz type of equation in E3, but the 
analysis also applies to the TE case. 

The approach in this paper leads to the spectral analysis of a quadratic opera- 
tor pencil, where the spectral parameter relates to the quasimomentum and where 
the frequency is a parameter. We show that the underlying operator has a discrete 
spectrum, where the eigenvalues are symmetrically placed with respect to the real and 



imaginary axis (Hamiltonian symmetry). A finite element method is used to discretize 
the operator pencil and the resulting matrix pencil is transformed into a gyroscopic 
pencil. We implement the structure preserving skew-Hamiltonian isotropic implic- 
itly restarted Arnoldi algorithm SHIRA [3TJ to compute eigenvalues of the gyroscopic 
matrix pencil. 

2. Bloch solutions and spectral problems. Let F denote the lattice 2tt1? 
and e a measurable function with 

e(x + 'y,u>) = e(x,u), V 7 £ T (2.1) 

for almost all x — (x\,Xi) 6 R 2 . A function with the property (|2.1[) is called Y- 
periodic. The fundamental periodicity domain of the lattice V is Vt = (—ir,ir) 2 and 
= 47r 2 denotes the area of the domain. In physics, a fundamental domain is often 
called a Wigner-Seitz cell [TBI [2]. We assume that e is in L°°(R 2 ) for each oj £ [ojq, us{\ 
and that e is invertible, positive and bounded 

< c < e(x,u) < ci, (2.2) 

for almost all x £ R 2 and w £ [wq,u>i]. Let T 2 = R 2 /T denote the torus in two 
dimensions and C°°(T 2 ) the space of smooth complex- valued T-periodic functions. 
The space L 2 (T 2 ) is defined to be the completion of C°°(T 2 ) in the L 2 -norm and 
H 1 ^ 2 ) is the completion of C°°(T 2 ) in the iJ^norm. 

Let V denote the gradient with respect to x £ R 2 and let V- denote the divergence. 
The transverse magnetic polarized waves (0, 0, E3, Hi, i?2, 0) can be written in terms 
of the electric field E alone. The equation 

-Ao = w 2 £(i,u)» (2.3) 

models a time-harmonic monochromatic electromagnetic wave with frequency u> and 
electric polarization E(x) = (0,0, v(x)). A Bloch solution of (|2.3p is a non-zero 
solution of the form 

v(x) =e ik - x u(x), (2.4) 

where u is a L-periodic function and k £ C 2 is the quasimomentum (or the Floquet- 
Bloch wave vector) [T7]. Since V(e lk ' x u(x)) = e lk x (\7 + \k)u(x) the Bloch solutions 
of (|2.3|) are T-periodic solutions of 

- (V + ifc) • (V + ik)u(x) = lu 2 e(x, w)u. (2.5) 

The frequency w as a multi-valued function of the quasimomentum k is called the 
dispersion relation and the graph of the dispersion relation defines the Bloch variety 

Definition 2.1. The complex Bloch variety is the set 

B = {(k,u>) £ C 2 x C \the problem h2.3\) has a non-zero Bloch solution} (2.6) 

Analytic properties of the Bloch variety B are discussed in Kuchment [T7]. The 
Bloch solutions of (|2.3p can be determined from spectral problems. In the frequency 
independent case e = e(x), Bloch solutions with k £ R 2 are given by eigenvectors 
of a selfadjoint operator. In the frequency dependent case e = e(x,to), the spectral 
problem is non-linear, which complicates the analysis. We will consider solutions 
(k,oj) £ C 2 x R in the frequency dependent case and use the real frequency u as a 
parameter. The two cases are discussed below. 
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2.1. The frequency independent case. Assume that e = e(x) independent of 
the frequency to and define the operator 



A(x,V)v = — ^-Av. (2.7) 

e{x) 

The spectrum of elliptic operators with periodic coefficient has been studied in detail 
by many authors including Odeh and Keller [53], Reed and Simon [21], and Figotin 
and Kuchment |6J. The L2-spectrum a(A) has a band structure and the existence 
of band-gaps for certain geometrical structures has been proved [6l [7] [8] . From the 
Floquet-Bloch theory [2UH7] follows the spectrum of A in the infinite periodic medium 
can be obtained from the k - dependent family of T-periodic solutions (|2.5p . We give 
below the well-known construction of the band structure. By definition, the function 
u £ i? 1 (T 2 ) is a solution of the shifted problem 

A(k)u = uj 2 u, A(k)u = — ^(V + ifc) • (V + ifcWaO (2.8) 

e(x) 

if the integral identity 

/ (V + ik)u • (V + i/c)0dx = uj 2 / u^edx (2.9) 
Jt 2 Jt 2 

is satisfied for all <j> 6 i? 1 (T 2 ). Each operator A(k) is selfadjoint in L 2 (T 2 ), with scalar 
product weighted by e(x) and k € M 2 . Since A(k) is non- negative and has a compact 
resolvent, the spectrum of A(k) consists of eigenvalues 

< w 2 (fc) < Lo\{k) < ... 

where the spectral parameter w 2 represent the square of the time frequency u of the 
wave. The main spectral result is that the spectrum of A (|2.7p is the union 

a(A)=\JS n , S n = [rmnu 2 n (k),maxLJ 2 (k)}. (2.10) 

n>0 

The intervals S n are called the stability zones of the operator A and the complement 
of cr(A) is called the instability zone. The successive segments S n and S n +i may 
overlap, but we have a gap in the spectrum cr(A) if the intersection S n n 
empty for some n. That is, if uj 2 belongs to the instability zone for all k 6 ffi. 2 . The 
spectrum cr(A) can also be described in terms of the Bloch variety (|2.6p . The real 
Bloch variety is the intersection 

B,cai = Bni 3 (2.11) 

and the spectrum of A is the projection of the real Bloch variety onto the u>- axis. 
This alternative description of the spectrum is important in the frequency dependent 
case below. 

3. The frequency dependent case. We assumed above that the permittivity 
e is independent of the frequency uj. The spectrum <r(A) coincides in the frequency 
independent case with the range of the dispersion relation uj 2 (k). The permittivity can 
be constant in a frequency interval, but only vacuum is independent of the frequency. 
If we exclude vacuum, the Kramers-Kronig relations [3, 13 implies that the imaginary 
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part of e(u>) cannot be zero for all to. We focus on the dispersive case e = e(x,u>), 
where e is real-valued in the frequency interval [woj^i]- That is, we assume that 
there is no absorption line in the frequency interval under consideration. This is a 
reasonable model for many dialectic materials (insulators) for a range of frequencies in 
the micro- wave regime and in the optical regime [T3j. Below, we study the quadratic 
pencil in the case of a real-valued permittivity function. 

3.1. The quadratic operator pencil. We presented in Section [2~T1 results for 
the shifted operator (|2.5p in the case e = e(x), where the spectral parameter uj 2 
corresponds to the square of the time frequency uj. We consider below the frequency 
dependent case and relate the spectrum parameter to the quasimomcntum k. We 
reduce the spectral problem in k to a problem in the complex amplitude of the vector 
k. 

Let k = Ak, where A £ C and A; is a unit vector in R 2 . Define a family of spectral 
problems as 

- (V + iXk) • (V + i\k)u(x) = oj 2 e(x, u)u, (3.1) 

where u and e are T— periodic functions. This spectral problem in (A, u) corresponds 
in the frequency independent case to the problem A(k)u — uj 2 u in (|2.8|) . but we 
consider k(oj) and not uj(k). The family of spectral problems (|3.1|) is quadratic in the 
spectral parameter A, with parameter uj. We write below the classical problem (|3.1| 
in a weak form. Define the sesquilinear forms 

a Q : i/ 1 (T 2 ) x i/ 1 (T 2 ) ->-C 

ai : L 2 (T 2 ) x H 1 ^ 2 ) -^C (3.2) 

a 2 : L 2 (T 2 ) x L 2 (T 2 ) -> C 

where 



ao(u, v) = / Vit • Vw — u 2 euv dx, 
Jt 2 

ai(u,w) = 2i / uk-Vvdx, (3-3) 
Jt 2 

a,2(u, v) = uvdx. 
Jt 2 

The weak solutions of p.ip are defined as eigenvectors u of i/ 1 (T 2 )\{0} and eigen- 
values A G C which satisfy 

A 2 02(m, v) + Xai(u, v) + ao{u, v) = 0, (3-4) 

for all v e H l (f 2 ). 

Lemma 3.1. The sesquilinear forms oq, a\ and a 2 are bounded for all elements 

u,v e H 1 ^ 2 ). 

Proof. The sesquilinear forms ao and a 2 are trivially bounded and the bounded- 
ness of ai follows from the inequality \ab\ < 6\a\ 2 /2 + \b\ 2 /(25), S > 0. □ 

Let (-, -) H i(T2) denote the scalar product in i? 1 (T 2 ). The spectral problem (|3.4p 
can be written as 



\ 2 (A 2 u, v) H i(j2) + \{Aiu, v) H i(j2) + (A u,v) H i^2 ) =0, (3.5) 
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where the operators Aq, A\ and A 2 are defined by the bounded sesquilinear forms 

(A fe u,«) H i( X 2) = a k (u,v), k = 0,1,2. (3.6) 

According to the Riesz theorem these operators are on iJ x (T 2 ) unique, linear and 
bounded [15]. The equality (|3.5j) that holds for all v e iJ x (T 2 ) is equivalent to the 
operator equation 

A Q u + XA lU + X 2 A 2 u = (3.7) 

in ff^T 2 ). 

Let £(_ff 1 (T 2 )) denote the set of all bounded linear operators on i? 1 (T 2 ). We 
introduce the operator pencil Q : C — > £(7? 1 (T 2 )) 

Q(A) = A a + XA X + A 2 A 2 , AeC. (3.8) 

The quadratic eigenvalue problem is then: Find A £ C and a non-zero u e i/ 1 (T 2 ) 
such that 

Q(A)u = 0. (3.9) 

The adjoint operator pencil is 

Q*(\) = A* + \Al + \ 2 A* 2 , AeC (3.10) 

and the pencil Q(X) is said to be selfadjoint if Ao, A\ and A 2 are selfadjoint [2D]. The 
lemma below shows that the pencil Q(X) is selfadjoint. 

Lemma 3.2. The operators Aq, Ai and A 2 are self-adjoint 

Proof. Since the operators Aq, Ai and A 2 are defined on the full space if 1 (T 2 ) 
the lemma follows from the equalities 

(A u,v) H i {T 2 } = a (u,v) = (A v,u) Hl( j2^ = {v, A^u) H1{T2) = (A^u, w) H i (T 2) 



(Aiu,v) H i(T2) = ai(u,v) = (A 1 v,u) H1{T 2 ) = (v, AJu) ffl(T2) = (AJu,u)hi(t=), 



(A 2 u,w) ff i (T 2 ) = a 2 (u,v) = (A 2 v,u) Hl(j2) = {v,A* 2 u) H1(J2) = {Alu,v) H x {r iy 

□ 

3.2. Spectral properties. We consider in this section spectral properties of 
the operator pencil (|3.8p . Fredholm theory is used to prove that the spectrum of Q 
is discrete, where the eigenvalues are symmetrically placed with respect to the real 
and imaginary axis. Moreover, we prove that the pencil Q(X) cannot be reduced to a 
monic bounded pencil, that is, to a pencil of the form 

Q(A) = B + \B 1 + A 2 I, AeC. (3.11) 

Nonmonic pencils has in general a more complex spectral structure |20j . 
Lemma 3.3. The operator A 2 has an unbounded inverse. 
Proof. Since A 2 > we have [2J 



^2w||jji(T2) < \\A 2 \\c(H 1 (T 2 ))(A 2 U,u) H i^2- ) 



= C / |w| 2 dx. 
Jt 2 
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(3.12) 



The operator Ai has a bounded inverse if and only if there exist a positive constant c 
such that ||A2u||hi(t 2 ) > c||-u||hi(t 2 ) for all u £ i? 1 (T 2 ). The inequality (|3.12|) implies 
that such constant c cannot exist (take for instance u(x) = e m ' x , n G N), but A2 has 
an inverse since the null space 7V(^2) is trivial. □ 

A bounded operator T is Fredholm if the range R(T) is closed and if the null space 
N(T), and the complement of the range R(T) are finite dimensional. The Fredholm 
index of a Fredholm operator T is defined as 

indT = nulT-defT, (3.13) 

where the nullity nulT is the dimension of N(T) and the deficiency defT is the 
dimension of the complement to range R(T). The resolvent set p is the set of all 
A G C, such that the operator T(A) is continuously invertible and the spectrum is 
defined by the complement 

a(T)=C\p(T). (3.14) 

The essential spectrum cr css (T) of T is the set of all A G cr(T) such that T(A) is not a 
Fredholm operator. The complement of the essential spectrum 

(T) = a{T)\a ess {T) (3.15) 

is called the discrete spectrum. The discrete spectrum consists of eigenvalues of finite 
geometrical multiplicity, nulT < 00. We shall show that the spectrum of Q{\) is 
discrete. The proof consists of two steps. We prove that the operator pencil Q{\) is 
a sum of compact operators and a Fredholm operator. 

Lemma 3.4. The operators A\ and A2 are compact in H l (Y 2 ). 

Proof. We show that A\ is compact. Let {u n } C if x (T 2 ) be a given bounded 
sequence. Then since if 1 (T 2 ) is compactly embedded into L 2 (T 2 ) 29J, the sequence 
{u n } has a convergent subsequence {u rn } in i 2 (T 2 ). Let u m and u m i denote two 
element in the subsequence {u m }. From the boundedness of the Hermitian form a\ 
follows 

\\Aiu m - Mu m '\\ 2 H nj2) — \ai{u m - u m i , A x u m — Aiu m i)\ 

yo. ID j 

< C\\Axu m - Aitt m '|| H i( T 2)||'U rn - u m /|| £ 2 (T 2). 

Then since {u m } is a Cauchy sequence in L 2 (T 2 ), the sequence {Aiu m } is a Cauchy 
sequence in iJ x (T 2 ), which converges since iJ 1 (T 2 ) is complete. Using the inequality 
||u||l 2 (t 2 ) ^ C|| u ll_ff 1 (T 2 )) t ne compactness of A2 is proved in the same way. □ 

Theorem 3.5. The operator pencil Q(X) is a Fredholm operator with index zero. 



Proof. The operator Aq in (|3.8p can be written of the form 

(A u, v) h i {1 2 ) = (A^u, v) H i(j2) + (A[, 2) u,w) h i (T 2), (3.17) 
where A^ and A^ are defined by the bounded sesquilinear forms 

(A^u, i>)iri(T 2 ) = — lo 2 / euwdx, (A^u, v)hi-(t 2 ) = / Vu • Vwdx. (3.18) 

JT 2 JT 2 

The operator A^ is selfadjoint and compact since e G L°°(T 2 ) is real- valued and 
i/ 1 (T 2 ) is compactly embedded into L 2 (T 2 ). Represent u, v G i/ 1 (T 2 ) with its Fourier 
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scries: 



i{x) = u n j n ' x , v(x) = J2 v m e im - x . (3.19) 



Since the gradients of u and v are in L (T ) the Plancherel's identity [32] gives 
(A^u, v) H i {T 2 ) = / V int^e" 1 '* • V 



= |«| ^ \n\ 2 u n v:, 



(3.20) 



n€Z 2 



which implies that the dimension of the null space is one. The dimension of the 
orthogonal complement of the range is also one, since i?(r)^ = N(T) for self-adjoint 
operators. That is, A Q is a Fredholm operator with index (13.13|) zero. The self- 
adjoint operator pencil Q(X) is a sum of A and compact operators, which is a 
Fredholm operator with the index of A^ [H] . □ 

3.2.1. Hamiltonian symmetry. We show in this section that the eigenvalues A 
of (|3.4[) are symmetrically placed with respect to the real and imaginary axes. Lemma 
13.21 implies that the spectrum of Q(X) is symmetric with respect to R, since Q(X) is 
invertible if and only if Q*(X) is invertible. 

Lemma 3.6. The operators Aq, A\ and A 2 have the properties Aqu — Aqu, 



Aiu = — A\u and A 2 u = A 2 u. 

Proof. The lemma follows from the equalities 



(A u,v) H i(T2) = {A u,v) HH jv ) = a (u,v) = (A u,v) H i(j2), 



(Aiu,v) H i( T 2) = {A 1 u,v) H1(J2) = -ai(u,v) = -(A 1 u,v) H i {1 2 ) , 



(A 2 u,v) H i {T 2 ) = (A 2 u,v) Hl(j2) = a 2 (u,v) = (A 2 u,v) H i {T 2 ) . 

□ 

The theorem below states the Hamiltonian structure of the spectrum a(Q(X)). 

Theorem 3.7. The spectrum of the quadratic operator pencil Q(X) consists of 
eigenvalues of finite multiplicity that have the Hamiltonian structure (A, —A, A,— A). 

Proof. Theorem ()3.5p states that Q(X) is a Fredholm operator, which implies that 
the spectrum o~(Q(X)) consists of eigenvalues of finite multiplicity. Lemma 13.61 gives 
the relation 



Q(X)u = A u + XAm + X A 2 u = A u - XAxll + X A 2 u = Q(~X)u, (3.21) 
which implies that —A is an eigenvalue whenever A is an eigenvalue. The two relations 
Q*(-A) =A -XA 1 + X 2 A 2 = Q(-A), (3.22) 

Q*{\) = A + XA 1 +X 2 A 2 =Q(A), (3.23) 
follows from Lemma |3~21 □ 



3.3. Band gaps. Assume that the periodic permitivity e is smooth in R 2 . If the 
equation 



Lv = 0, L = A + u 2 e(x,u) (3.24) 

has a bounded solution \v(x)\ < C, then (|3.24p also has a Bloch solution with a real 
quasimomentum k. The existence of a Bloch solution with k G M 2 implies that the 
operator L is not invcrtiblc in L2(R 2 ), see Chapter 4 in [T7]. We will use the above 
results to define a band-gap in terms of the Bloch variety. 

Assume that the parameter uj in Q(X), that corresponds to the time frequency, is 
real. The complex Bloch variety of the pencil Q(A) is the set 

B{Q) = {(k,uj) G C 2 x R\Q(X)u = Ohas a non-zero solution}, (3.25) 

and the real Bloch variety of the pencil Q(X) is the intersection 

Brcai(Q)=5(Q)ni 3 . (3.26) 

We state below the condition for a band-gap in terms of the real Bloch variety. 

Definition 3.8. The frequency lu is in a band gap if the projection of the real 
Bloch variety onto the to-axis is empty. 

Solutions of the form (k,ui) eI 2 xR exist if the frequency ui is not in a band 
gap. That is, bounded Bloch waves exist for the given frequency u and we say that 
u> belongs to the stability zone. A gap in the spectrum of Q(X) is also called a 
polarization gap since we only consider the polarization E(x) — (0,0, u(x)). Notice 
that we let k be real in the frequency independent case and defined a band gap in 
(|2.10|) as an empty intersection of two successive bands 

[mino; 2 (fc),maxw 2 (fc)] n [mincj 2 +1 (fc), maxw, 2 1+1 (fc)] = {}. (3.27) 

We also used that the spectrum alternatively can be described in terms of the projec- 
tion of the real Bloch variety (|2.11[) onto the w-axis. 
The dual (or reciprocal) lattice to T = 2nJ? is 

r* = {q G M 2 | 7 - q G 27rZ,V 7 G T} (3.28) 

We define the fundamental domain (the Brillouin zone) of the dual lattice T* — 1? as 
the set 



2' 2 



(3.29) 



The real part of all eigenvalues A to (|3.ip is CI*- periodic, since the Bloch solutions 
p.4p are periodic in the real part of k — kX. We write the complex spectral parameter 
A in the form 

A = A r + Aii, (3.30) 
where A r and Aj are real numbers. The set 

Bn* = {(A r fc, Aifc, w)el]*xl 2 xM Q(X)u = has a non-zero solution}, (3.31) 



is a subset of the Bloch variety (|3.25|) . This set corresponds to solutions with the real 
part of the wave vector k = kX in the Brillouin zone f2* and a real time frequency to. 
The frequency u> is in a band-gap if Aj for all eigenvalues of Q(X), with fcA r € CI*. 

When u> = and k = (1,0), the lemma below shows that the eigenvalues can be 
calculated explicitly. 

Lemma 3.9. Let oj = and k = (1,0). Then all eigenvalues A £ C can be written 
of the form 

A = Z + iZ. (3.32) 

Proof. Let v(x) = v m e lm ' x , to £ Z 2 and represent u 6 iJ 1 (T 2 ) with its Fourier 
series 

u{x) = {L r l e m ' x - (3.33) 

n<£Z 2 

The equation (Q(A)m, u)#i(T2) = gives 

(m + Afc) • (m + \k)u m i>* m = 0. (3.34) 
The coefficients are zero whenever u rn ^ which implies 

f Iml 2 + 2A r fc • m + A 2 - A 2 = 0, 

< A ' (3.35) 

Ajfc • to + A r Ai = 0. 

Let A; = (1,0) and assume A; = 0. Then follows A = A r = —mi and to,2 = 0. The 
secound equation in (|3.35[) is satisfied for A r = —mi and Ai = ±|m-2| follows from the 
first equation in (|3.35p . □ 

Notice that eigenvalues A in the lemma above is real or purely imaginary if to-2 = 0. 
This case corresponds to a permittivity e that only dependence on one coordinate. 
We illustrate the lemma in Section [5] 

The frequency u) — is an eigenvalue of Aq = Aq(lo) and the corresponding eigen- 
vector is constant. An eigenvalue to of the operator- valued function Aq(lh) corresponds 
to a solution of the classical eigenvalue problem 

— Ait = u) 2 e(x, uj)u, (3.36) 

where u and e are T-periodic. Assume that {u, Co) is a solution of Aq{uj)u = 0. The 
pencil Q(A) can then be written of the form 

Q(X)u = X(A lU + XA 2 u). (3.37) 

Zero is always an eigenvalue of the pencil and a nonzero eigenvalue A is a solution of 

A lU + XA 2 u = 0, (3.38) 

where u e N(Aq(uj)). All eigenvalues of (|3.37p are real and the number of eigenvalues 
at u> = lu depends on the dimension of the null space N(Aq(uj)). 

4. Numerical treatment. 
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4.1. Discretizing the problem. In this section we discuss the numerical com- 
putation of approximate eigenvalues of the operator pencil Q{X). According to The- 
orem [23 the spectrum of Q(X) has a Hamiltonian structure. A numerical method 
should preserve the symmetry of the spectrum. We will show that this can be achieved, 
if we discretize the weak formulation of the problem, which was given in (|3.4p , by a 
conforming Galerkin ansatz. 

Let Vjv be an iV-dimensional linear subspace of iJ 1 (T 2 ), and let {</>i, . . . , tp^} 
be a basis of Vjy. The subspace Vjv can be constructed, for example, using a finite 
element ansatz. We seek Xn G C and um S Vjv, such that 

X 2 N a 2 (u N ,v N ) + X N ai(u N ,v N ) +a (u N ,v N ) = 0, (4.1) 

for all vn £ V/v- The approximate eigenfunctions are of the form 

N 

u N = ^2a j <f> j , (4.2) 

where the coefficients ctj are given by the solution of the quadratic matrix eigenvalue 
problem 

A^A 2 u + AjvAiu + A u = 0, (4.3) 
with u = (cti, . . • , ccn) T ■ The N x N matrices A 2 , Ai and Ao are given by 

(A ) m „ = / V</>„ ■ V(j) m - io 2 e(uj)(j) n (j) m dx, 

(Ai) mn = 2i / <j) n k ■ V(f> m dx, (4.4) 

JT 2 



(A 2 ) mn = / 0„0m dx, m,n = l,...,N. 

JT 2 

Since Vat is a subspace of 7f 1 (T 2 ), we observe that the matrices in (|4.3p have the same 
properties as the operators in the operator pencil 



A 2 =A 2 H >0, A 1= A", A =A£. (4.5) 

Note that the matrices A 2 and Ao are real, while the matrix Ai is purely imaginary. 
Since we want to avoid complex matrix arithmetic in actual computations, we trans- 
form the problem as follows: Letting Xn = — i/ij K = — Ao, G = iAi and M = A 2 
the quadratic matrix eigenvalue problem can be rewritten as 

/i 2 Mu + /iGu + Ku = 0, (4.6) 

The N x N matrices in (|4.6|) are real and satisfy 

M = M T > 0, G = -G T , K = K T . (4.7) 



The quadratic matrix eigenvalue problem in (|4.6p is therefore called a gyroscopic ma- 
trix eigenvalue problem. This name refers to the fact, that quadratic matrix eigenvalue 
problems of this form typically occur in the analysis of undamped gyroscopic systems 
[3Tj . It is easy to show, that the set of eigenvalues of a gyroscopic eigenvalue problem 
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has a Hamiltonian structure: Suppose that fj, is an eigenvalue of (|4.6p . Then, there 
exists an eigenvector u £ C N , u ^ 0, such that 

/i 2 Mu + /iGu + Ku = 0. (4.8) 

Due to the structure of the matrices in (|4.6p . we have 

7Z 2 Mu + JZGu + Ku = |i 2 Mu + /iGu + Ku = 0, 
(-n) 2 Mu - /xGu + Ku = (/j 2 Mu - itGu + Ku) T = 0, (4.9) 
(-Jt) 2 Mu - JZGu + Ku = (/i 2 Mu - /iGu + Ku) H = 0. 

From ()4.9|) follows that every element of the set {fi, ~p, — //, — jtl} is an eigenvalue. Hence, 
the discrete problem will preserve the Hamiltonian structure of the continuous prob- 
lem. This statement remains true as long as the discrete space Vm is a subspace of 
iJ^T 2 ). In the context of finite clement methods this means that the finite element 
basis functions satisfy periodic boundary conditions. 

4.2. Linearization. The numerical solution of quadratic and especially of gy- 
roscopic matrix eigenvalue problems is still an active field of research P32 [SSJ IS] • A 
comprehensive overview of numerical algorithms that can be applied to quadratic 
matrix eigenvalue problems can be found in [31] . 

A common way to treat quadratic matrix eigenvalue problems is to transform 
the quadratic problem to an equivalent generalized eigenvalue problem of double di- 
mension 9, 3 lj . The gyroscopic matrix eigenvalue problem in (|4.6[) can be linearized 
as 

Ax = /iBx, (4.10) 
where the 27V x 2N block matrices A and B are given by 

A-(-K -o). *=(l m). *=(,!)■ <"» 

Since M is positive definite and I is the identity matrix, the block matrix B is also 
positive definite. The computation of the eigenvalues /x can be done, for example, 
with the QZ-algorithm [3"T] . 

It turns out, however, that this algorithm is not well-suited for the compuation 
of photonic band gaps, mainly for two reasons: First, the matrices are typcially large 
and sparse. This is due to the fact, that they are generated by a finite element 
ansatz. The QZ-algorithm will destroy the sparsity in general and therefore lead to 
excessive storage requirements. Second, only a few eigenvalues need to be computed in 
order to determine the photonic band structure. For these two reasons Krylov space 
methods, such as an implicitely restarted Arnoldi process (IRA) or an implicitcly 
restarted Lanczos process (IRL), are preferable Using shift-and-invert strategies 
these algorithms can be used to compute a small number of eigenvalues fi\,. ..,Hm 
that are closest to a prescribed shift parameter fj,Q. Each computed eigenvalue \n 
determines up to four eigenvalues of the gyroscopic matrix eigenvalue problem, namely 
{He,~p e , -fit, -~Pi}- 

4.3. The SHIRA algorithm. A drawback of standard Krylov space methods 
is that they do not take into account the Hamiltonian structure of the set of eigenval- 
ues of the gyroscopic matrix eigenvalue problem. Suppose that /i £ C is an eigenvalue 

11 



for (|4.6p that is close to the real axis. A Krylov space method may find eigenvalues 
Hi,... , /xm, amoung which are approximations to ^ and ^ to /j, and to ~p, respec- 
tively. Numerical round-off errors could however lead to \it ^ JLgi . Thus, one would 
get a wrong extra quadruple of eigenvalues. 

As a remedy, one might consider to use structure-preserving algorithms [12; , like 
the SHIRA algorithm |22j . SHIRA stands for skew-Hamiltonian implicitely restarted 
Arnoldi process. As the name suggests, it is a modified Arnoldi process that preserves 
the Hamiltonian symmetry of the eigenvalues of a skew-Hamiltonian matrix. The 
algorithm can be used to compute a small number of eigenvalues for a large and 
sparse gyroscopic eigenvalue problem. Below we summarize the SHIRA algorithm 
and refer to [22] for details. 

The the gyroscopic matrix eigenvalue problem in (14.61) is linearized as 



Hy = M Sy, (4.12) 
where the 27V x 27V block matrices H and S are given by 

H = (m ~o K ) • s = ("o m) ■ y = (T ) 

The matrix H turns out to be a Hamiltonian matrix, i.e., it satisfies (HJ) T = HJ, 
where the 27V x 27V block matrix J is the imaginary unit matrix 

' - (_? 1) ■ 

The block matrix S is skew-Hamiltonian, i.e., it satisfies (SJ) T = — S3. The factor- 
ization 



s = v o 2 m){o \ ) ( 4 - 15 ) 

allows us to rewrite (|4.12p as 

Wy = /iy (4.16) 
where W is the 27V x 27V block matrix 

w-(S LTU -o k )C <«') 

The matrix W turns out to be Hamiltonian [22] . Since the spectrum of such a matrix 
has a Hamiltonian symmetry, the eigenvalues of the matrix 

R= (W-/i )- 1 (W + Ai0 )- 1 (W-^)- 1 (W + I^)- 1 (4.18) 

will have the same symmetry for every shift parameter /io £ C. Moreover, the eigen- 
values of R that are largest in modulus correspond to eigenvalues of W that are 
closest to elements of the set {/zo, — A*o,7*0: ~ Mo}- Therefore, an implicitely restarted 
Arnoldi process will in general converge to these eigenvalues. It should be noted that 
the application of R can be realized efficiently in actual computer programs. 

It can be shown that the block matrix R is skew-Hamiltonian. In exact arith- 
metic, the Krylov spaces 

Kn = span{v ,Rv ,R 2 v , . . . ,R n_1 vo} (4.19) 
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that are generated by R should therefore be isotropic, i.e., they should satisfy JC n _L 
J/C„. Due to numerical round-off errors, however, this property is usually lost af- 
ter some Krylov space iterations, and this will eventually destroy the Hamiltonian 
symmetry of the computed eigenvalues. In order to cope with this problem, SHIRA 
introduces an additional re-orthogonolization step, which ensures that the Krylov 
spaces remain isotropic up to machine precision. Because of this, SHIRA always con- 
verges to complex conjugate pairs of eigenvalues in the right hand complex halfplanc. 
The algorithm thus prevents the occurrence of wrong extra eigenvalue quadruples. 

The SHIRA algorithm has already been successfully used to compute the corner 
singularities in elasticity problems pQ. In our computations we found that in general 
it needs less implicit restarts than a standard implicitely restarted Arnoldi process. 
However, in some cases, especially when the relevant eigenvalues where close to each 
other, we could not obtain convergence with the SHIRA algorithm. In these cases we 
resorted to the linearization in (|4.10[) and used an IRA algorithm as provided by the 
ARPACK package. 

5. Numerical examples. To compute the band gaps, we choose frequencies ui 
in a frequency range [uj a ,ujb\ and let 

One can show that the band structure can be calculated reducing the angles to < 
9 < 7r/4, which is called the irreducible Brillouin zone [16 . The operator pencil Q(A) 
is discretized with quadratic elements, which results in a quadratic matrix eigenvalue 
problem on the form (|4.6p . Solving these quadratic matrix eigenvalue problems, we 
obtain for each frequency ui and each angle 9 a set of eigenvalues A(oj, 9). Eigenvalues 
A G A(cj, 9) that satisfy |A| < 1/ cos#, correspond to quasi-momentum vectors k = Xk 
that lie in the first Brillouin zone. A frequency u> is a band gap frequency if the 
imaginary parts of all eigenvalues in A(w, 9) do not vanisch for all 9. 

We compute band gaps when the periodic structure is a rectangular array of 
dielectric cylinders in air. The relative diameter of the cylinders with respect to the 
lattice constant is 0.75. Below we consider one frequency independent model and 
one frequency dependent model. The mesh width is h = 0.025 in both cases, which 
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Fig. 5.1. Imaginary part of A as a function of the frequency uj when k = (1,0). Left: Material 
parameters ei = 1 and t2 = 8.9. Right: t\ = 1 and e-x according to the model 



h/2n 


DOFs 


First gap 


Second gap 


Third gap 


0.050 
0.040 
0.030 
0.025 


2608 
4124 
7412 
10872 


(.24719, .27015) 
(.24710, .27016) 
(.24718, .27015) 
(.24711, .27015) 


(.41064, .45632) 
(.41064, .45631) 
(.41063, .45631) 
(.41029, .45650) 


(.61757, .66173) 
(.61755, .66171) 
(.61754, .66171) 
(.61754, .66171) 



Table 5.1 

Computed band gaps for the frequency independent model (Dobson). 



leads to quadratic matrix eigenvalue problems with approximately 10° unknowns. 
The method was implemented in Matlab and the computations were performed on a 
standard laptop. 

5.1. A frequency independent model. As a first illustration, we assume that 
the permittivity of the cylinders is 62 = 8.9, and that the permittivity of air is ei = 1 
for all frequencies under consideration. This model was numerically investigated in 
[3J. Since the permittivity is frequency independent, the band gaps can be calculated 
from the shifted problem (|2.8p . Dobson [3] used this formulation and calculated u>{k) 
where k is on the boundary of the irreducible Brillouin zone; See Figure I3~T1 It is know 
that a band gap generically opens on the boundary, but it is possible to construct a 
material function that gives a band gap opening inside the Brillouin zone [lOj . 

With the approach in this paper, a band gap corresponds to a frequency uj, where 
the imaginary part of all eigenvalues A are non-zero. Figure ISTTI shows the imaginary 
part of the eigenvalues X(uj, k) with smallest imaginary part when k = (1, 0) and 
uj G [0,0.7]. Notice that 9A is close to ±1 for uj close to zero. We know from Lemma 
13.91 that all solutions have SA <E Z when uj = 0. To compute band gaps, k is varied 
over the irreducible Brillouin zone with small steps in 9 and uj G [0,0.7]. The band 
gap frequencies in Table 15.11 match well with Dobson's results [3] . 

5.2. A frequency dependent model. We assume that the matrix material 
(air) has the permittivity t\ — 1 for all frequencies under consideration. The permit- 
tivity of the dielectric cylinders is generally a complicated function of the frequency. It 
is common in practice to use least-squares minimization to determine a rational func- 
tion approximation from frequency-domain data. The choice of the material model is 
not important for the algorithm, since the frequency dependence is only a parameter. 
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h/2n 


DOFs 


First gap 


Second gap 


Third gap 


0.050 
0.040 
0.030 
0.025 


2608 
4124 
7412 
10872 


(.28231, .30031) 
(.28230, .30031) 
(.28231, .30031) 
(.28228, .30031) 


(.44251, .47728) 
(.44251, .47728) 
(.44250, .47727) 
(.44238, .47784) 


(.60209, .63320) 
(.60208, .63318) 
(.60207, .63318) 
(.60207, .63318) 



Table 5.2 

Computed band gaps for the frequency dependent model (5. 2(1 . 




e/n 



Fig. 5.2. The first band gap tube 3A(lj, 8) for the material model (5.21) . 

We will consider the simple model 

e 2 (w)=a+— (5.2) 

where a = 1, b = 5.34, and c — 1. The constants in (|5.2p are chosen such that the mean 
value (e(0) +e(0.7))/2 is approximately 8.9. We assume that this model is valid for all 
frequencies in the range u> £ [0, 0.7]. This assumption make it possible to compare the 
frequency dependent model with the calculations for the frequency independent case 
above. Figure UTTl shows the imaginary part of the eigenvalues X(u>, k) with smallest 
imaginary part when k = (1, 0) and to s [0, 0.7]. A frequency u; is in a band gap when 
all eigenvalues have a non-zero imaginary part. The lowest band gap in Table 15.21 
corresponds in Figure f5T2l to the tube. The band gap can alternatively be illustrated 
with spectral surfaces w(6,X), where A is a real eigenvalue. The spectral surfaces 
in Figure 15.31 are computed with the presented method by sorting the computed 
eigenvalues X(ui, k). 

6. Discussion and conclusions. We have analyzed a method to compute Bloch 
waves in periodic and frequency dependent materials. The problem is formulated as a 
quadratic operator pencil in the quasimomentum k. This formulation is, in the case of 
frequency dependent materials, a significant simplification of usual approach, where 
the time frequency u> is calculated as a function of k. The pencil can be linearized since 
the non-linearity is polynomial. This is important, because many robust numerical 
algorithms for computing eigenvalues are only available for linear eigenvalue problems. 
The problem was discretized with finite elements and approximate eigenvalues was 
successfully computed with an implicitly restarted Arnoldi process. 
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Fig. 5.3. The spectral surfaces 10(8, A) for the material model (5.21) 
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